#'###########################################
#'                                          #
#'  ------------------------------------    #
#'  |     COIN     __                  |    #
#'  |            >(' )                 |    #
#'  |              )/                  |    #
#'  |             /  `----/            |    #
#'  |             \  ~=- /             |    #
#'  |   ~^~^~^~^~^~^~^~^~^~^~^~^~^~^   |    #
#'  ------------------------------------    #
#'                                          #
#'###########################################


rm(list=ls())

#'#
#'# Chemins de travail
#'#
R_ROOT <- "C:/UQAR/Recherche/Maitrise/Donnees"
R_WORKING_DIRECTORY <- file.path(R_ROOT, "R")


#'#
#'# Chargement des dependances
#'# 
source(file.path(R_WORKING_DIRECTORY, "Utils/loadDependencies.R"))
DATA_PATH <- file.path(R_RESULTS_PATH, "HRCompare/results")

hr.stats.raw <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_raw.csv")), type="raw", stringsAsFactors=FALSE)
hr.stats.LC3 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_LC3.csv")), type="LC3", stringsAsFactors=FALSE)
hr.stats.LC32 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_LC32.csv")), type="LC32", stringsAsFactors=FALSE)
hr.stats.LC321 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_LC321.csv")), type="LC321", stringsAsFactors=FALSE)
hr.stats.mb1 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_mb.csv")), type="mb1", stringsAsFactors=FALSE)
hr.stats.mb2 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_mb2.csv")), type="mb2", stringsAsFactors=FALSE)
hr.stats.radius250 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_radius250.csv")), type="radius250", stringsAsFactors=FALSE)
hr.stats.radius5000 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_radius5000.csv")), type="radius5000", stringsAsFactors=FALSE)
hr.stats.radius5000_n50 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_radius5000_n50.csv")), type="radius5000_n50", stringsAsFactors=FALSE)
hr.stats.filter <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_filter.csv")), type="filter", stringsAsFactors=FALSE)
hr.stats.rect250x1000 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_rect250x1000.csv")), type="rect250x1000", stringsAsFactors=FALSE)
hr.stats.rect8000 <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_rect8000.csv")), type="rect8000", stringsAsFactors=FALSE)
hr.stats.tracksGen <- cbind(read.csv(file.path(DATA_PATH, "HR_compare_tracksGen.csv")), type="tracksGen", stringsAsFactors=FALSE)


hr.stats <- rbindlist(list(hr.stats.raw, hr.stats.LC3, hr.stats.LC32, hr.stats.LC321,
				hr.stats.mb1, hr.stats.mb2, hr.stats.radius250, hr.stats.radius5000, hr.stats.radius5000_n50, hr.stats.filter, hr.stats.rect250x1000,
        hr.stats.rect8000, hr.stats.tracksGen))

hr.stats[, type := factor(type, levels=c("raw", "LC3", "LC32", "LC321",
            "mb1", "mb2", "filter", "tracksGen", "radius250",
            "rect250x1000", "radius5000", "radius5000_n50",
            "rect8000"))]

setkey(hr.stats, type, method, percent)

hr.stats <- hr.stats[sizeRatio < 200]

statsCols <- c("sizeRatio", "overlapIndex", "propInRef")
stats <- hr.stats[, c(lapply(.SD, mean, na.rm=TRUE), lapply(.SD, sd, na.rm=TRUE)), by=list(type, method, percent),
    .SDcols=statsCols]
setnames(stats, 4:9, c(paste(statsCols, "mean", sep='.'), paste(statsCols, "sd", sep='.')))
setkey(stats, method, percent)


hr95 <- stats[J(c("mcp", "kernel"), c(100, 95))]

hr95 <- hr95[!type %in% c("LC3", "radius5000_n50")]

cols <- grey(c(0.8, 0.4))

g1 <- barchart(overlapIndex.mean ~ type, data=hr95, groups=method, col=cols,
         scales=list(rot=45), origin=0,
         ly=hr95[, overlapIndex.mean - overlapIndex.sd],
         uy=hr95[, overlapIndex.mean + overlapIndex.sd],
         prepanel=prepanel.errorbar, panel=panel.errorbar,
          key=list(text=list(c("Kernel",
                     "MCP")), 
             rectangles=list(col=cols), x=1, y=0.98, corner=c(1, 1))
         )

g2 <- barchart(sizeRatio.mean ~ type, data=hr95, groups=method, col=cols,
   scales=list(rot=45), origin=0,
   ly=hr95[, sizeRatio.mean - sizeRatio.sd],
   uy=hr95[, sizeRatio.mean + sizeRatio.sd],
   prepanel=prepanel.errorbar, panel=panel.errorbar,
   key=list(text=list(c("Kernel",
               "MCP")), 
       rectangles=list(col=cols), x=1, y=0.98, corner=c(1, 1)))


g3 <- barchart(propInRef.mean ~ type, data=hr95, groups=method, col=cols,
    scales=list(rot=45), origin=0,
    ly=hr95[, propInRef.mean - propInRef.sd],
    uy=hr95[, propInRef.mean + propInRef.sd],
    prepanel=prepanel.errorbar, panel=panel.errorbar,
    key=list(text=list(c("Kernel",
                "MCP")), 
        rectangles=list(col=cols), x=1, y=0.98, corner=c(1, 1)))

print.multi(list(g1, g2, g3), hfill=T,
    save=TRUE, save.dest="C:/UQAR/Recherche/Maitrise/Donnees/Resultats/compareHRs.png")


filt <- hr95[type %in% c(c("LC32", "radius5000", "radius5000_n50", "rect8000", "tracksGen"))]
g4 <- barchart(sizeRatio.mean ~ type, data=filt, groups=method, col=cols,
 scales=list(rot=45), origin=0,
 ly=filt[, sizeRatio.mean - sizeRatio.sd],
 uy=filt[, sizeRatio.mean + sizeRatio.sd],
 prepanel=prepanel.errorbar, panel=panel.errorbar,
 key=list(text=list(c("Kernel",
             "MCP")), 
     rectangles=list(col=cols), x=1, y=0.98, corner=c(1, 1)))
     
print.multi(list(g4), hfill=T,
    save=TRUE, save.dest="C:/UQAR/Recherche/Maitrise/Donnees/Resultats/compareHRs2.png")

filt <- hr.stats[type %in% c(c("LC32", "tracksGen"))]
filtStats <- filt[, c(lapply(.SD, mean, na.rm=TRUE), lapply(.SD, sd, na.rm=TRUE)), by=list(type, method, percent),
    .SDcols=statsCols]
setnames(filtStats, 4:9, c(paste(statsCols, "mean", sep='.'), paste(statsCols, "sd", sep='.')))
setkey(filtStats, method, percent)

filtStats[, method2 := paste(method, percent, sep="_")]
barchart(sizeRatio.mean ~ type, data=filtStats, groups=method2, col=c("red", "blue"),
    scales=list(rot=45), origin=0,
    ly=filtStats[, sizeRatio.mean - sizeRatio.sd],
    uy=filtStats[, sizeRatio.mean + sizeRatio.sd],
    prepanel=prepanel.errorbar, panel=panel.errorbar)

ano <- aov(sizeRatio.mean ~ type + method2, data=filtStats)
summary(ano)


filt2 <- filt[, c(lapply(.SD, mean, na.rm=TRUE), lapply(.SD, sd, na.rm=TRUE)), by=list(type, method, percent, deployment),
    .SDcols=statsCols]
setnames(filt2, 5:10, c(paste(statsCols, "mean", sep='.'), paste(statsCols, "sd", sep='.')))
setkey(filt2, method, percent)
filt2[, method2 := paste(method, percent, sep="_")]
ano <- aov(sizeRatio.mean ~ type + method2 + type*method, data=filt2)
summary(ano)
plot(ano$residuals)
#'#
#'# Anova et Test de Tukey
#'# 
ano <- aov(overlapIndex.mcp.mean ~ type, data=hr.stats)
summary(ano)
tk <- TukeyHSD(ano)
plot(tk)

#'#
#'# Analyse de puissance
#'# 
f <- effect.size.anova(hr.stats, overlapIndex.mcp.mean, quote(type))
p <- pwr.anova.test(k=7, n=20, f=f, sig.level=0.05)
